1. Field of the Invention
The invention relates generally to the filtered backprojection technique for 3D image reconstruction in a cone beam x-ray imaging system, and more specifically to the division of the acquired cone beam detector data into groups for individualized filter processing.
2. Description of the Prior Art
A filtered backprojection (FBP) cone beam image reconstruction technique is described by Kudo, H. and Saito, T., in their article entitled "Derivation and Implementation of a Cone-Beam Reconstruction Algorithm for Nonplanar Orbits", IEEE Trans.Med, Imag., MI-13 (1994) 196-211.
Briefly, the FBP technique consists of the following steps at each cone beam view (i.e., at each position of the radiation source as it scans about the object, and at which an imaging detector acquires a corresponding set of projection data):
1. Compute a ID projection (i.e., line integral) of the measured cone beam image acquired on a detector plane 1 at each of a plurality of angles .theta.. This step is illustrated by FIG. 1A for a given angle .theta..sub.1 of a plurality of angles .theta., where the projection 2 at coordinates (r,.theta.) comprises the integrated values of the cone beam image 4 on detector plane 1 along a plurality of parallel lines L(r,.theta.) that are normal to angle .theta., each line L being at an incremental distance r from an origin O. PA1 2. Filter each 1D projection in accordance with a d/dr filter, resulting in a new set of values at each of the r,.theta. coordinates, such as shown by filtered projection 6 for the angle .theta..sub.1 in FIG. 1A. PA1 3. Normalize the filtered projections with a normalization function M(r,.theta.). Normalization is needed to take into account the number of times the plane of integration Q(r,.theta.) which intersects the source position and the line L(r,.theta.), intersects the scan path, since the data developed at each scan path intersection creates a contribution to the image reconstruction on the plane Q(r,.theta.). PA1 4. Backproject the filtered projection 6 from each angle .theta. into a 2D object space 7 which coincides with the detector plane 1. This step is illustrated by FIG. 1B, where in lines 8 spread the value from each r,.theta. coordinate into 2D space 7 in a direction normal to each .theta.. PA1 5. Perform a ID d/dt filtering of the backprojection image formed in 2D space 7 by step 4. The ID filtering is performed in the direction of the scan path, i.e., along lines 10, where the arrowhead points in the direction of the scan path. PA1 6. Perform a weighted 3D backprojection of the resulting data in 2D space 7 (i.e., from each pixel in the detector) onto a plurality of sample points P in a 3D object volume 12. The density assigned to each point P is weighted by the inverse of the square of the distance between the point and the spatial coordinates of the x-ray source (see Equation (59) of the forenoted Kudo et al article). PA1 1) Apply a mask to the set of cone beam projection data acquired by the detector at each of the source positions, so that only specific non-overlapping contributions to the Radon data can be developed from the projection data. PA1 2) Calculate line integral derivatives in the masked data. PA1 3) Perform a 2D backprojection of the derivative data onto an extended height virtual detector. PA1 4) Perform a 3D backprojection of the 2D data from the virtual detector into a 3D object space.
Generally, if the detector plane 1 comprises an N by N array of pixels, then the number of angles .theta. is typically given by .pi.N/2.
The above prior art procedure will be referred to hereinafter as the 6-step process. It is assumed in this process that the entire cone beam image of the object is captured on the detector of the imaging system. Consider a plane Q(r,.theta.), which intersects the object, formed by the source and the line L(r,.theta.) on the detector at angle .theta. and at a distance r from the origin. Ignoring the function M(r,.theta.), the operations 1 through 6 compute the contribution to the reconstructed object density on the plane Q(r,.theta.) from the x-ray data illuminating the plane and its immediate vicinity. Since the 6-step process is detector driven, a contribution from the data illuminating the plane is computed every time the plane intersects the scan path and thus is illuminated by the x-ray beam. Consequently, the function M(r,.theta.) is used after the filter function in step 2 to normalize the results. Normalization is particularly undesirable since it requires pre-computing and storing a 2D array M(r,.theta.) for each source position along an imaging scan path. Since there are usually hundreds, if not thousands of source positions, this type of normalization is both computationally intensive and resource (computer memory) expensive.
In U.S. patent application Ser. No. 09/052,281 entitled EXACT REGION OF INTEREST CONE BEAM IMAGING USING 3D BACKPROJECTION, filed Mar. 31, 1998 and incorporated herein by reference, inventor K. Tam departs from the conventional Radon space driven conversion processing techniques for image reconstruction (such as known by his U.S. Pat. Nos. 5,257,183 and 5,453,666), and discloses a way to incorporate the technique of data combination for region-of-interest (ROI) reconstruction, with the Kudo et al. image reconstruction processing, thereby providing an image reconstruction technique for a cone beam imaging system that can not only have a spiral scan path, but can also use a short detector, In K. Tam's technique, instead of division by the function M(r,.theta.) as done by Kudo et al., the effect of the normalization of the reconstructed object densities is achieved by dividing the x-ray beam coverage of integration plane Q(r,.theta.) between the various source positions that illuminate the plane without any overlap.
K. Tam's technique comprises a 4 step process:
The presence of a detector mask ensures that the contributions developed by processing projection data of the different detectors are unique and non-redundant (FIG. 1 of this disclosure). Accordingly, division by the function M(r,.theta.), or its equivalent, is no longer needed; a significant simplification in the image reconstruction signal processing. However, although step 2, is not complex, it is computationally expensive. More specifically, it comprises calculating a plurality of line integrals L(r,.theta.) on each set of the masked detector data, to generate sampled ID projections of the detector data. Line integral derivatives are then computed from the ID projections by taking the difference between parallel line segments L.sub.1 and L.sub.2, as shown in mask 200 of FIG. 2 herein. Note that the L.sub.1 and L.sub.2 line segments are not limited by the boundaries of the mask, and therefore their use results in an exact calculation for the derivatives of line integrals L(r,.theta.). This type of masking is referred to herein as "soft masking". Additional details of such soft masking can be found in K. Tam's recently issued U.S. Pat. No. 5,748,697, incorporated herein by reference. Step 3 backprojects the line integral derivatives onto the extended "virtual" detector. Before the 3D backprojection in step 4, the gradient of the backprojected virtual detector data in the direction of the scan path is calculated, and the result is then backprojected into the 3D object space for reconstruction the ROI of the object. For good image quality, the sampling of the projections and the number of source positions needs to be very fine.
Thus, the filter process described by this U.S. Ser. No. 09/052,281 is computationally costly.
In U.S. patent application Ser. No. 09/106,537 entitled SIMPLIFIED CONE BEAM IMAGE RECONSTRUCTION USING 3D BACKPROJECTION, filed Jun. 29, 1998 and incorporated herein by reference, inventor K. Tam introduces a Feldkamp convolution processing simplification (also referred to as ramp filtering) into the above-described image reconstruction processing, wherein the entire filter process of U.S. Ser. No. 09/052,281 is replaced with a simple single step of ramp filtering of the detector data in the direction of the scan path. This simplification is illustrated in FIG. 3, where L, L.sub.1 ' and L.sub.2 ' are three closely spaced parallel line segments that are bound by a mask 300, and L is midway between L.sub.1 ' and L.sub.2 '. Line segment L is representative of many such line segments formed at various angles in mask 300, and corresponds to the previously described lines L (r,.theta.) of FIG. 1, which as well known to those skilled in this technology are used for computing Radon derivative data from the cone beam projection data. In the technique described in U.S. Ser. No. 09/106,537, due to the bounding of the line segments L.sub.1 ' and L.sub.2 ' by mask 300, the Feldkamp convolution processing simplification (referred to as ramp filtering) is performed as a substitute for the line integral derivative calculations, which filter processing corresponds to calculation of the Radon derivative of the partial plane defined by the line segment L and the current source position, up to a multiplicative constant.
Although this operation is computationally very fast, it yields only an approximation of the Radon derivative of the partial plane, due to errors that come about due to the "hard masking" of the endpoints of line segments L.sub.1 ' and L.sub.2 ' by mask 300, as compared to the "soft" masking shown in FIG. 2. That is, it incorrectly limits the detector pixel values to those pixels that are in the mask area, and zeros out the detector pixel values that are outside of the mask boundaries, instead of correctly limiting only the line segments L to the mask area (and calculating the line integral derivatives using the unmasked original detector data when appropriate, i.e., near the mask boundaries).
Accordingly, in U.S. patent application Ser. No. 09/123,574 filed Jul. 27, 1998 entitled MASK BOUNDARY CORRECTRION IN A CONE BEAM IMAGING SYSTEM USING SIMPLIFIED FILTERD BACKPROJECTION IMAGE RECONSTRUCTION, K. Tam describes a technique for computing 2D correction data which, when combined with the ramp filtered 2D approximation data sets, yields an exact image reconstruction. As described in greater detail in U.S. Ser. No. 09/123,574, the mathematical difference between hard and soft masking, which involves only detector data around the mask boundaries, is calculated to arrive at an additive correction term. Although this technique is mathematically correct, artifacts are still possible due to a mismatches in spatial resolution and/or position between the correction data and the ramp filtered data. Furthermore, there is no apparent way to eliminate the artifacts from this "approximation+correction" technique.
It would be desirable to provide an image reconstruction technique that has a processing speed closer to that of the forenoted U.S. Ser. No. 09/106,537, but with the image reconstruction accuracy of the forenoted U.S. Ser. No. 09/052,281.